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ABSTRACT: A very fruitful approach to the solution of image segmentation 
and surface reconstruction tasks is their formulation as estimation problems 
via the use of Markov random field models and Bayes theory. However, the 
Maximum a Posteriori (MAP) estimate, which is the one most frequently 
used, is suboptimal in these cases. We show that for segmentation problems, 
the optimal Bayesian estimator is the maximizer of the posterior marginals, 
while for reconstruction tasks, the thresholded posterior mean has the 
best possible performance. We present efficient distributed algorithms for 
approximating these estimates in the general case. Based on these results, 
we develop a maximum likelihood procedure that leads to a parameter-free 
distributed algorithm for restoring piecewise constant images. To illustrate 
these ideas, the reconstruction of binary patterns is discussed in detail. 
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1. Introduction. 


A very powerful approach for solving computational early vision problems is 
to formulate them as estimation problems. The estimates are based on probabilistic 
models that embody both our prior generic knowledge about the behavior of 
the solution, as well as the process by which the observations are obtained. In 
particular, various researchers have used Bayesian estimation and Markov Random 
Field (MRF) models to perform image segmentation (Elliot et. al. (1983); Hansen 
and Elliot (1982); Geman and Geman (1983); Marroquin (1984b); Cohen and 
Cooper (1984)) and surface reconstruction tasks (Marroquin (1984a)). Most of these 
researchers have assumed that the Maximum a Posteriori (MAP) estimate, which 
is defined as the configuration which maximizes the posterior distribution, provides 
the best possible solution to these problems, and several methods (such as dynamic 
programming and "simulated annealing" (Kirkpatrick, et. al. (1983)), a form of 
stochastic relaxation) have been proposed to compute it. 

Let us consider, however, the example portrayed in figure 1. 

Panel (c) represents the MAP estimate of the binary MRF (a) (whose parameters 
are assumed to be known) from the noisy observations (b). It is clear that the 
estimates shown in (d) and (e) (which we will describe later) are better than the MAP 
estimate from almost any viewpoint. This suggests the use of criteria other than 
the value of the posterior probability to measure the performance of an estimator. 
In particular, we propose the use of the following criteria for the two classes of 
problems mentioned above: 

1.1. Error Criterion for the Segmentation Problem. 


Consider a field x with N elements each of which can belong to one of a finite 
set Qi of classes. Let x { denote the class to which the i ih element belongs. The 
segmentation problem is to estimate x from a set of observations {y 1} ..., y p }. Note 
that xi does not necessarily correspond to the image intensity. It may represent, for 
example, the texture class for a region in the image (as in Elliot et. al., 1983), etc. 

A reasonable criterion for the performance of an estimate x is the number of 
elements that are not classified correctly. Therefore, we define the segmentation 
error e a as: 

N 

e s (x } x) = (1 — 6(x{ — 2 t -)) , X{, Xi G Qi (1) 

*=i 


where 


f 1, if a = 0 

10, otherwise 
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(d) (e) 


Figure 1. (a) Sample function of a binary MRF. (b) Output of a binary symmetric channel 
(error rate: 0.4) (c) MAP estimate, (d) Monte Carlo approximation to the MPM estimate, (e) 
Deterministic approximation to the MPM estimate. 

1.2. Error Criterion for the Reconstruction Problem. 

In this case, we also consider a field x with N elements which can take values 
on finite sets {<5,}, but now we assume specifically that x ,• represents the intensity 
of an image (or the height of a surface) at site i. This suggests that an estimate x 
should be considered "good” if it is close to x in the ordinary sense, so that the 
total squared error: 

MM) = ]C(z *-*») 2 ® 

t=i 

will be a reasonable measure for its performance. 

Let us now derive the optimal estimators for these error criteria. 
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2. Optimal Bayesian Estimators. 


In this section we introduce two estimators based on the posterior distribution 
P x \ y : The Maximizer of the posterior marginals x MVhU which is obtained by 
maximizing separately the posterior marginal distribution of each element, and the 
Thresholded posterior mean x TPM . Their formal definition is as follows: 


= q : Pi\y(q; y) = sup {P t | y (r; y)} (4) 

r£Q; 


where P^ y is the posterior marginal distribution of the i ih element: 

Pi\y{qj y) — X^ Px\yi x \y) (5) 

x:xi=q 


The TPM estimate is: 


x T PM{i) = q : (9 - Zi') 2 = {(r - z t ) 2 } (6) 


where x is the posterior mean: 


5 = y&v) (7) 

x 

We will show that the MPM and TPM estimators are optimal with respect to the 
error criteria (1) and (3) respectively. In what follows, we will assume that the prior 
and conditional distributions P x , P y \ x are known, so that the posterior distribution 
can be obtained from Bayes rule. 

Proposition 1: x MPM is optimal with respect to the segmentation error e s . 

This is not a new result (see for example Abend, 1968); however, we present here 
a simple proof: 

Let x = xmpm, and let 2 denote any other estimate for x. We will show that 

E[e s [x,x)] < E[e s (x,z)] 

where the expectation is taken over all x and all possible observations y. 

For any fixed y we have, from (1): 

N 

e s (x | y ) = E[e s {x, x) | y] = £ £ (i - <5(z t - - 5 t -))P, |y (x; y) 

X t'=l 

N N 

X/ X/ Px\y{ x > y) ~ XT(^ ~ Pi\y{ x ii 2/)) 

'_ 1 _ / * - ’ 
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with P^y given by (5). Also, 


N 


M z \ v) = EO - p i\y( z i\y)) 

i— 1 


Therefore, 


N 


e,(i | y) - e s (z | y) = - E [ Wl!/) - p ,^ v)] 

i= 1 



But by the definition (4) of x, each term on the right hand side of (8) is non-positive, 
so that 

e s (x | y) - e s (z | y) < 0 for all y 


and therefore, 


E [ e s[x,z)] < E[e a (x,z)) | 


Note that we have not made any assumption about P x , P y]x , or {Q t }, so that 
this is a completely general result. In particular, the sets of valid classes Q { do not 
have to be the same for all i , so that this result holds if, as in Geman and Geman 
(1983), x is the pair {/,/}, where / denotes the intensity of a piecewise constant 
image, and l the boundaries between uniform regions. 

Proposition 2: x TPM is optimal with respect to the reconstruction error e r . 

Proof: Again, we put x = x TPM , and let z denote any other estimate. For any fixed 
y we have, from (3): 


where 


also, 


M* I y) = E E (*< ~ X if p x\,(x; y ) 

x 1 = 1 
N 

= IT Yl( X i “ 2x i x i + *i) P z\yfa y) 

i=l x 

= E [(*.' - *i? + - 5. ? ] 

1=1 

*i = E X i p x\,fa v) 

X 

= T, x ! P x\ v {x; y ) 

z 

M Z I y) = K 5 *' “ Z if + X \ - x i] 

i=l 
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so that 


N 

t r (x | y) - e T (z I y) — [(5, - z,-) 2 - [x x - z x f) < 0 

i=l 

from the definition (6) of x i 

Remarks: 

1. In the limit of fine quantization the variables x x become continuous (i.e., Q x =R 
for all i). In this case we recover the classical result that establishes that the posterior 
mean minimizes the mean squared error; for the particular case of Gaussian (Habibi, 
1972; Nahi and Assefi, 1972) or isotropic fields (Levy and Tsitsiklis, 1983), it is then 
possible to construct efficient algorithms that compute this quantity. 

2. The surface reconstruction problem treated by Marroquin (1984a) can be 
considered a particular instance of the reconstruction problem defined here. If we 
let x be the pair {/, /}, where / is the (continuous valued) depth process and l is 
a binary line process, the optimal estimator will be x TPM (note that for a binary 
variable, x TPM {i) = x MPM (i)). 

3. These results can also be applied to the case of an m-ary line process (like the 
one suggested by Geman and Geman (1983)) by defining an appropriate mixed 
error: 

e m(f, l, f, l) = (/* “ ft ) 2 + ^ TAl- &(lj — k)) 

t=l 1=1 

with S as defined in equation (2). One can show, using arguments similar to those 
of the proofs above, that the optimal estimate, for any positive value of X will 
^ {Itpm^mpm}- Note that the posterior mean for / and the posterior marginal 
distributions for l must be computed by averaging over all possible values of both 
/ and l: 

/ ‘ 

p iM) = E L Pf,i\ y {fJ\y) 

1 W,'=5 


3. Algorithms. 

As in the case of the MAP estimate, the exact computation of the optimal MPM 
and TPM estimates is a formidable task, except for very small N. It is possible, 
however, to define general distributed Monte Carlo procedures that will allow us to 
approximate these values as precisely as we want. 

The general idea is to use the Metropolis algorithm (Metropolis, 1953) to 
generate global configurations of the field 2 which are, asymptotically, distributed 
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according to P x | y. Since the Markov chain defined by the Metropolis algorithm is 
known to be ergodic (Geman and Geman, 1983), we can approximate x by: 


n — k 


n 

E 

t—k 


.(*) 


and the posterior marginals by: 



n 


k — n 


E - 9) 


t—k 




where x® is the configuration generated by the Metropolis algorithm at time t, and 
k is the time required for the system to be in thermal equilibrium. From these 
values, x M pm and xtpm can be easily computed using (4) and (6). 

If we compare this procedure with the use of’’Simulated annealing” (Kirkpatrick, 
et.al M 1983) for approximating the MAP estimate, we note that even from a 
computational viewpoint it has some distinct advantages: 

1. It is difficult to determine in general the descent rate of the temperature 
(annealing schedule) that will guarantee the convergence of the annealing 
process in a reasonable time (it usually involves a trial and error procedure). 
Since we are running the Metropolis algorithm at a fixed temperature (T = 1), 
this issue becomes irrelevant 

2. Since in our case we are using a Monte Carlo procedure to approximate 
the values of some integrals, we should expect a nice convergence behavior, in 
the sense that coarse approximations can be computed very rapidly, and then 
refined to an arbitrary precision. 

The main disadvantage of this procedure is that in the case of the segmentation 
problem, a large amount of memory might be required if the number of classes 
per element m is large (we need to store the N(m - 1) numbers that define the 
posterior marginals). 

With respect to the relative performance, we point out that for high signal to 
noise ratios, the MAP estimate is usually close to the optimal one (which explains its 
good performance in many cases). If the noise level is high, however, the difference 
in the performances of the two estimators may be dramatic (see fig. 1). 

Finally, from a conceptual viewpoint, the definition of the solution of the 
estimation problem in terms of the posterior mean (or the posterior marginals) 
may lead in some cases to the construction of efficient deterministic algorithms for 
approximating the solution. To make these ideas concrete, we now analyze in detail 
a specific instance of the segmentation problem. 


6 



4. Optimal Segmentation of Binary MRF’s. 


In this section we will consider the following particular example of the general 
segmentation problem: 

Let / represent a binary pattern on a rectangular lattice L with N elements, 
which is a particular realization of a Markov Random Field with potentials given 
by: 


r-1, if fi = fj and ||i-j|| 6(0,1* 

~ 1. if k^fj and ||i-ill e (0,1 

to, otherwise 

where i and j are sites of L and ||i - j'|| is the distance betweeen i and j (i.e., / 
is a realization of a two-dimensional lsing net). Without loss of generality, we will 
assume that /,• 6 {0,1} for all i e L. 

The probability density of the configurations F is given by: 

Pr(F = f)=l exp|-i-t/(/)] (11) 

where Z is a normalizing constant; the energy U(f) is; 

m= e nfi.fi) 

*» jENn 


where i,j £ N n means that i and j are nearest neighbors, and T 0 is the natural 
temperature of the system. 

Suppose that / is sent through a binary symmetric channel with error rate e, 
so that the conditional distribution of the observations g is: 


P(9i | /.) = {£ e) ' 

The posterior distribution is: 


if 9% = fi 

if 9i 7 A fi i /«» 9% € {0,1} 




( 12 ) 


where Zp is a constant, and the posterior energy Up is: 

u r = jrY, v (fi ,/;) + “E(1- S{fi - s,-)) (13) 

0 » 


with 

a = ln ( i 7 i } 
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1.17 



Figure 2. Ratio of the average errors of die MAP and MPM estimators for a 2 x 2 Ising net. 
4.1. Error Analysis. 

Using equations (1) and (12), it is not difficult to write down the exact expressions 
for the average segmentation error corresponding to the MAP and MPM estimators 
for this problem. However, the numerical evaluation of these expressions is feasible 
only for small values of N. 

In figure 2 we show a plot of the ratio: 

_ ^[ e Af/LPj 
E[ e AtPMl) 

for a 2 x 2 lattice, for different values of the error rate e and the natural temperature 
To. As expected, r is never less than 1. In the worst case (for e — 0.1 and To = 0.2) 
the error of the MAP estimate is 1.17 times that of the MPM estimate; if T 0 is not 
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too small and c is not too large, both estimates coincide, and as c approaches 0.5 
(low signal to noise ratio), the MPM estimate is consistently better than the MAP. 
An experimental analysis of larger lattices reveals a similar qualitative behavior, but 
the values of r are much larger in this case. 

4.2. Example. 

We now return to the example presented in figure 1, and examine it in more 
detail. Panel (a) represents a typical realization of a 64 x 64 lsing net with free 
boundaries, using a value of T 0 = 1.74 (0.75 times the critical temperature of the 
lattice); panel (b), the output of a binary symmetric channel with error rate e = 0.4; 
panel (c) the MAP estimate, and panel (d) an approximation to the MPM estimate 
(which we will label "MPM (M.C.)") obtained using the Metropolis algorithm and 
equation (10) to estimate the posterior density. The corresponding values of the 
posterior energy Up (equation (13)) and the relative segmentation error (e s /64 2 ) are 
shown on table 1. 


Table 1 


f 

9 

A 

/map 

fMPKfiM'C.) 

A 

Energy 

-5594.8 

-226.0 

-6660.9 

-6460.0 

-6427.0 

Seg. Error 

— 

0.4 

0.33 

0.128 

0.124 


4.3. A Fast Algorithm. * 

For this problem, it is possible to construct a fast deterministic algorithm whose 
experimental performance (in terms of the average segmentation error) is equivalent 
to the Monte Carlo method discussed above. It is based on the following ideas; 

First, we note that for a binary pattern, the MPM and TPM estimates coincide. 
We will approximate the posterior mean of (12) by that of a Gaussian distribution 
Pg with the property: 


Pa{h) = Pe-WV for all h e {0,1}. 

Lp 


In particular, we use: 


P G (h) = j- exp[-± £ (h { - htf - a UK ~ 9i ) 2 }. 

* G l0 hj£N N t 


For this distribution, h is the (unique) minimizer of the convex function: 

Uc(h) = ~ E (hi-htf-aUlH-K? 
l0 ijeN N i 
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which corresponds to the unique fixed point of the system: 



Eyeyv, hf 1 + aT 0 gi 
m+ecTo 


where 

N t = {jeL : ||x - ill = 1>. 
We could now approximate our estimate by putting: 


ft = e(S,-) 



where 

G(z) = I 1 ' lfl >5 

10, otherwise 



There is an additional consistency condition that / must satisfy, however. It can 
be shown that when the posterior distribution has the form given by (12) and (13), 
the MPM estimate /, which by definition satisfies: 




also satisfies: 

Pi\g{fi) f) > Pi\g({l ~ (16) 

which means that if we replace the observations by the MPM estimate, and compute 

a new MPM estimate for this modified problem, we should get the same result. 

* 

Translating this condition to the case of / , we get that it must satisfy: 


u = ©(O 

where h‘ satisfies: 

. _ Eye*,- h) + aToQ(h') 

|AT,-| + aT 0 

In practice, we get h' as the fixed point of the system: 

,-11,11 &€», ‘l 1 * 1 + oWlk' 1 * 1 ) 

*•-rarrsi- 




with 

h"® = 1 
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Note that the function: 


UkW = E (*.• - h jf + *t 0 Z(h, - ef^j) 2 

i,j€Ni i 


acts as a Lyapunov function for the system (18), which is therefore (locally) stable 
(Vidyasagar, 1978). 

This algorithm can be visualized as operating in two steps: In the first one, 
we extract all the information that we need from the observations and encode it in 
h (which is continuous-valued), and in the second one, we find the closest binary 
pattern that satisfies the consistency condition (16). 

A * 

To illustrate the performance of this approximation, we show f , for the 
example discussed above, in panel (e) of figure 1 , and its corresponding energy and 
segmentation error in the last column of table 1 (labeled "MPM det."). 

4.4. Dimensionality of the Parameter Space. 


Equations (14) and (18), which describe the deterministic approximations to 
Impm depend on the parameters of the system, £ and T 0 , only through the product: 


7 = aTo = To log 



which means that the behaviour of the algorithm is completely characterized by the 
single parameter 7 . In the case of the Monte Carlo approximation (10), if we fix 
the value of 7 , the value of To cannot be chosen arbitrarily, since it has to satisfy 
the consistency condition: 


with 



a . 1 1 ”£ 

or = log ( ~Y~ , , 


1 N 


where 2 is the residual process defined as: 


2 , = < 1, iffiJ^9i 
: 0 , otherwise 


( 20 ) 


( 21 ) 
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This means that, given 7 , the correct value of T 0 can, in principle, be determined 
in an adaptive way, so that in this case too, the behaviour of the approximation 
depends effectively only on 7 . 

5. Simultaneous Estimation of the Field and the Parameters. 

The above considerations suggest that we can estimate simultaneously the 
original field and the parameters of the system by looking for the supremum of an 
appropriate likelihood function over the one-dimensional subspace parametrized 
by 7 . We proceed as follows: 

For a given value of 7 , we can approximate the corresponding MPM estimate 
f using the methods developed in the previous section, and compute the residual 
process z and the conditional (on 7 ) Maximum Likelihood Estimate of the error 
rate e using equations (20) and (21). The corresponding conditional estimate for T 0 
will be: 

% = l ( 22 ) 

To measure the "likelihood” of the estimate /, we use the degree of uniformity 
(or "whiteness") of the residual process 2 . This property can be quantified by the 
variance of the local noise density, which we estimate as follows: 

We cover the lattice with a set {S 3 } of m non-overlapping squares (say, 8 
pixels wide). For each square Sj, the relative variance of the noise density is: 



with 




where \Sj\ is the area of the j th square. 
The desired likelihood function is: 




Note that a more conventional likelihood function, such as the conditional 
likelihood proposed by Besag (1972), will not work in this case; this function is 
defined as: 


L{f) = 


with 


LkO) = n p 0i \fjJeNi,f 0 ) = 

*ECk 
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_ exp l-jrLieN'VCfi'h)] _' 

ex p| — fo ^ eN - W + eXp i _ R V{1 - fi,fj)\ 

= IKi + explr- E V 0r,ljW\ * = 1,2 

i(zCk To 

where the "codings" C\ and C 2 are sets of conditionally independent sites defined 
as: 



C\ = {i : (i t - is odd and y{ is even) or (x t is even and y t - is odd)} 

C 2 — {i : [z{ is odd and t/ t is odd) or (z t - is even and y t - is even)} 

with ( xi,yi ) denoting the row and column indices of site i. In our case, we find 
that as 7 decreases, / becomes more and more uniform, while t 0 remains almost 
constant. It is not difficult to see that as a result, the conditional likelihood L will 
decrease monotonically with 7 , which renders it useless for our purpose. 

The range of values [70,7 m ] of the parameter 7 that corresponds to the class 
of systems of interest can be determined as follows: 

One can show that for 7 > 8 we will always have 9i for all h so that 

we can use 7 ^ = 8. The value of 70 can be obtained from an upper bound for e 
and a lower bound for T 0 . For example, assuming that 6 < .45 and T 0 > . 5 T c , we 
get 70 = .23. (Note that when the natural temperature T 0 of a first order, isotropic 
MRF is below 0.5 times T c (the critical temperature of the lattice; see Kindermann 
and Snell, 1980), the patterns become practically uniform (i.e., =constant for all 
2 ), while for values of To greater than 1.5T C , we get patterns that are practically 
indistinguishable from white noise. Therefore, the assumption T 0 > . 5 T C covers 
practically all the interesting cases). 

The complete estimation procedure is as follows: 

Maximum Likelihood Estimation Algorithm: 

1: Sample the interval [ 70 , 7 ^] at the points 


7o ^ 7i> • • *7 n ^ 7 M 


2: For each 7 £ © = {7i> • • -7n} : 

2.1: Find /( 7 ) using (14) and (18). 

2.2: Compute 2 using (21). 

2.3: Compute e using (20). If e = 0 , set I(/( 7 )) = — 00 and proceed with the 
next value of 7 . Otherwise, compute a and go to 2.4. 

2.4: Compute to using (22). 
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2.5: Compute using (23) and (24). 

A * 

3: Compute the optimal estimate f using: 

/ = f(l) : Uf{ 7*)) = sup £(/( 7 )) (25) 

7 eQ 

* * * 

The corresponding e , T 0 will be the optimal estimates for e and T 0 , respectively. 
Remarks: 

1. This estimation algorithm allows us to reconstruct a binary pattern / from 
the noisy observations g without having to adjust any free parameters. The only 
prior assumptions correspond to the qualitative structure of the field / (first order, 
isotropic MRF) and to the nature of the noise process, but neither the natural 
temperature T 0 nor the error rate e have to be known in advance. In practice, this 
means that we can apply it to restore any binary image with uniform granularity, 
even if it has not been generated by a Markov random process. We have used this 
algorithm to reconstruct a variety of binary images with excellent results. In figure 3 
we show such a restoration. The observations (b) were generated from the synthetic 
image (a) with an actual error rate of .35 (assumed unknown). The MLE for / is 
shown in (c). 

2. This procedure can be easily extended to handle any one-parameter noise 
corruption process (such as zero mean, additive white Gaussian noise). The extension 
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to the case of N-ary fields, i.e., to the restoration of piecewise constant images, is 
also straightforward (using the general algorithm (10) together with equation (4) 
instead of (14) and (18) in step 2.1). 

3. We have found that the likelihood function (24) is reasonably well behaved as 
a function of 7 . This permits us to perform the one-dimensional search for its 
supremum in an economical way, by first determining its approximate location using 
a coarse sampling pattern, and then refining its position by a fine sampling of a 
reduced interval. In practice, it is possible to get very good results using on the 
order of 15 samples. 

4. The whole procedure is highly distributed, so that it is possible to implement it 
in parallel hardware in a very efficient way. 

6. Conclusions. 

A very fruitful approach to the solution of image segmentation and surface 
reconstruction tasks is their formulation as estimation problems, via the use of 
MRF models and Bayes theory. However, the MAP estimate, which has been often 
used, will not always give the best solution. We have shown that, for segmentation 
problems, the optimal Bayesian estimator is the maximizer of the posterior marginals 
(MPM), while for reconstruction tasks, the thresholded posterior mean (TPM) has 
the best possible performance. In some cases, these estimators will give results 
that are similar to those obtained by the MAP procedure, but in many interesting 
situations (in particular for low signal to noise ratios) their performance will be 
significantly better. 

The exact computation of these estimates (including the MAP) is computationally 
unfeasible in most practical situations. We have presented a general Monte Carlo 
procedure which approximates their value as precisely as desired (at the expense of 
increased computational cost). This procedure can be easily implemented in parallel 
hardware, and it has the advantage that coarse solutions are computed very fast, 
and are then progressively refined. We have also shown how in particular cases, it is 
possible to use the conceptual framework provided by these estimators to construct 
fast algorithms with very good experimental performance. 

Finally, we have developed a maximum likelihood procedure for the simul¬ 
taneous estimation of a piecewise constant field and the parameters of the system. 
This construction leads to a parameter-free distributed algorithm for reconstructing 
piecewise constant images from noisy observations. 
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